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Abstract. The formation and evolution of binaries which 
contain two neutron stars or a neutron star with a black 
hole are discussed in detail. The evolution of the distribu- 
tions in orbital period and eccentricity for neutron star bi- 
naries are studied as a function of time. In the model which 
hts the observations of high mass binary pulsars best the 
deposition of orbital energy into common envelopes has to 
be very efficient and a kick velocity distribution has to con- 
tain a significant contribution of low velocity kicks. The 
estimated age of the population has to be between sev- 
eral 100 Myr and 1 Gyr. The birthrate of binary neutron 
stars is ~ 3.4 x 10~ 5 yr~ 1 (assuming 100% binarity) and 
their merger rate is ~ 2 x 10~ 5 yr -1 . The merger rate of 
neutron star binaries is consistent with the estimated rate 
of 7-ray bursts, if the latter are beamed into an opening 
angle of a few degrees. We argue that PSR B2303+46 is 
possibly formed in a scenario in which the common enve- 
lope is avoided while for the other three known high-mass 
binary pulsars a common envelope is required to explain 
their orbital period. 



Key words: binaries: close - stars: neutrons - pulsars: 
general - Galaxy: stellar content 

1. Introduction 

Peculiarly enough, the first binary pulsar to be discov- 
ered, PSR B1913+16 (Hulse & Taylor 1975), turned out 
to contain two neutron stars one of which is visible as a 
recycled pulsar and the other is not visible. The majority 
of about thirty recycled pulsars discovered so far in the 
Galactic disc are found to have white dwarf companions. 
Four of these systems may have a neutron star companion 
instead, as is suggested by their mass. 

An adequate model of the population of Galactic neu- 
tron star binaries is highly desirable, as their mergers 
are currently considered as the most viable sources of 
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detectable bursts of gravitational wave radiation (Clark 
& Eardley 1977) and possibly 7-rays (Paczyhski 1986). 
Estimates of the merger rates of the Galactic binaries 
which contain two neutron stars based on observations 
are strongly biased by the uncertainties in the pulsar 
lifetime, beaming factor and the fraction of unobserved 
low-luminosity pulsars (Curran & Lorimer 1995; van den 
Heuvel & Lorimer 1996). Therefore, it is useful to inves- 
tigate the Galactic population of neutron star binaries by 
simulating their formation and the subsequent orbital de- 
cay due to gravitational wave radiation and to fit the or- 
bital parameters to observations. The population of binary 
neutron stars is a part of the more general family of de- 
scendants of massive binaries. Hence, an adequate model 
has to comply with observations also in the stages preced- 
ing the formation of neutron stars (see Meurs & van den 
Heuvel 1989 for an example of this approach). 

After the pioneering studies of Clark et al. 1979, Ko- 
rnilov & Lipunov (1984a, b) and Dewey & Cordes (1987) 
the formation of single and binary neutron stars and the 
merger rate of the latter are modelled by, among oth- 
ers Tutukov & Yungelson (1993a, b), Portegies Zwart & 
Spreeuw (1996), Lipunov et al. (1996, 1997), Bagot (1996). 

Since the paper by Flannery & van den Heuvel (1975), 
a lot of effort was spent on reproduction of the parame- 
ters of PSR B1913+16. But surprisingly, the rest of the 
family didn't attract much attention. The primary aim of 
most studies was to constrain the orbital parameters of im- 
mediate progenitors of high-mass binary pulsars (HMBP, 
hereafter) and kicks imparted to neutron stars (see e.g. 
Wijers et al. 1992, Yamaoka et al. 1993; Fryer & Kalogera 
1997; Bagot 1997; Lipunov et al. 1997). No attempts were 
made to restrict the models by simultaneous confrontation 
with all observed systems and to understand how typical 
they are in the total inferred population of double neutron 
stars. Though only four HMBP are definitely known 1 a 

1 We omit the globular cluster object B2127+11C which 
could be formed by a dynamical encounter rather than by 
binary evolution (Phinney & Sigurdsson 1991) and B1820-11 
which may be young or have a non-neutron star as companion 
(Phinney & Verbunt 1991). 
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better understanding of their formation is helpful, as it 
may provide additional constraints on the probable rate 
of events accompanying their merger. 

In Sect. 2 we discuss the formation scenarios for 
HMBP. (We use the terms "neutron star" as a generic 
term for the remnant of a high mass star and "pulsar" 
for a neutron star which is observable in radio.) Partic- 
ular attention is given to the revision of the "standard" 
scenario suggested by Brown (1995). Section 3 gives the re- 
sults of our population synthesis computations and the in- 
ferred orbital-period-eccentricity diagram for binary neu- 
tron stars. This probability diagram provides a unique op- 
portunity to compare the combined results of evolutionary 
computations and population synthesis with observations 
of HMBP. A requirement for pulsars to be located in the 
regions of this diagram with high likelihood also allows to 
constrain the progenitor population and its evolution. We 
derive the estimate of the merger rate of the Galactic neu- 
tron stars from the model which fits best to the observed 
population. Discussion follows in Sect. 4 and conclusions 
are summarized in Sect. 5. 

2. Formation of binary neutron stars 

2.1. Scenarios of formation 

Neutron stars are believed to be born in the supernova 
explosions triggered by the iron-core collapse in massive 
stars. The range of progenitor masses of neutron stars in 
isolation assumed in the present paper is 8 to 40 M . How- 
ever, the lower limit may be actually even higher than 10 
M©; Ritossa et al. (1996) demonstrated that a 10 M Q star 
ends its evolution as 1.26 M© oxygen-neon white dwarf. 
This would require a reduction of the model occurrence 
rate of supernovae from single stars by ~ 30%. If a star 
is stripped from its hydrogen envelope due to Roche-lobe 
overflow in a close binary the lower mass limit for the 
formation of a neutron star increases to about 12 Mq. 
The necessary condition for the formation of a neutron 
star is the requirement for the helium star remnant of 
a binary component to have a mass in excess of about 
2.2 to 2.5 M Q (Tutukov & Yungelson 1973; Habets 1985), 
which enables the formation of a carbon-oxygen core of a 
Chandrasckhar mass. The upper limit is less well known. 
The main factors contributing to this are uncertain stellar 
mass-loss rates, poor understanding of mixing processes 
in stellar interiors, the processes during the supernova it- 
self and the uncertainty concerning the equation of state 
of the compact object which results from the supernova 
(see e.g. Brown et al. 1996). The "standard" assumption is 
that neutron stars are formed in binaries by stars initially 
less massive than 40 M© (van den Heuvel & Habets 1984). 
However, this limit may be as high as ~ 60 M© if helium 
stars have very strong stellar winds (Woosley et al. 1995). 
On the other hand, there are indications from theory as 
well as from observations that this limit may be as low as 



20-25 M Q (see Portegies Zwart et al. 1997b and Timmes 
et al. 1996) The wide spread in the mass of observed can- 
didate black holes in binaries suggests that there may be 
factors other than the progenitor mass alone which deter- 
mines the fate of a star (Ergma & van den Heuvel 1997). 

The chain of events resulting in the formation of a 
bound pair of neutron stars is known now for more than 
20 years, starting from studies by van den Heuvel & 
Heise (1972), Tutukov & Yungelson (1973), De Loore et 
al (1975), Flannery & van den Heuvel (1975) and van den 
Heuvel (1976). Let us introduce the notation ms, gs, he, 
co, ns and bh, for a main-sequence star, giant, stripped 
helium core or helium star, the stripped carbon-oxygen 
core of a helium star, neutron star and a black hole, re- 
spectively (see Portegies Zwart & Verbunt 1996 for an ex- 
tended description). A bullet above a stellar type indicates 
that this star is to experience a supernova explosion in the 
next evolutionary step. A subscript e means that the orbit 
of the system is eccentric. A bound pair is indicated by 
symbols enclosed by round brackets: two main-sequence 
stars in a circular orbit are denoted as (ms, ms) and an ec- 
centric binary neutron star as (ns, ns) e . A square bracket 
by the symbol indicates that this star fills its Roche lobe 
and stably transfers mass onto its companion. If the bi- 
nary experiences a phase of unstable mass transfer the 
symbols are placed within braces. Using this notation the 
"standard" scenario for the formation of an (ns, ns) e is: 

I. (ms, ms) e , [gs, ms), (he, ms), [he, ms), (co, ms), (ns, ms) e , 
{ns, gs}, (ns, he), {ns, he}, (ns, co), (ns, ns) e . 

Actually, all other scenarios for the formation of HMBP 
differ from scenario I in the number of Roche-lobe over- 
flow and common envelope events experienced by the bi- 
nary. The stages [he, ms) and {ns, he} typically are absent 
for helium stars more massive than ~ 4 M Q which do not 
expand much beyond ~ 10 Rq after core helium exhaus- 
tion (Paczyhski 1971; Habets 1985; Woosley et al. 1995). 
The strong stellar wind of helium stars may also prevent 
the occurrence of both phases of mass transfer from the 
helium star. 

A peculiar case are binaries with a primary star ini- 
tially more massive than <~ 40 M© and a secondary less 
massive than 25 to 40 M©. They may avoid the first 
Roche-lobe overflow at all due to their strong stellar wind 
which may 1) prevent them from expanding to a radius 
exceeding several hundred solar radii and blow away most 
of their envelopes during the Wolf-Rayet phase, and 2) 
increase the orbital separation via loss of angular momen- 
tum. However, the secondary star may overflow its Roche 
lobe when it evolves onto the giant branch. The orbital 
period subsequently reduces as a result of the spiral-in 
and the second supernova occurs in a short-period binary. 
Note that a kick can also decrease the orbital separation 
upon the formation of the first neutron star. 
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Scenario I involves stages when the neutron star is im- 
mersed into a common envelope (the {ns, gs} stages and, 
on some occasions, also {ns, he}). Recently certain doubts 
were expressed concerning this phase based on the work 
of Zel'dovich et al. (1972, see also Chevalier 1993; Brown 
1995; Fryer et al. 1996) which have shown that neutrino 
cooling allows hypercritical accretion onto neutron stars, 
which is O(10 8 ) higher than the canonical photon Edding- 
ton limit of ~ 10~ 8 M yr _1 . Common envelopes in mas- 
sive binary systems may provide suitable conditions for 
hypercritical accretion; the net outcome may be the col- 
lapse of the neutron star into a black hole. Whether or 
not the binary survives the common envelope and a (bh, 
ns) e binary forms or the spiral-in continues until both 
stars merge into a single object is nor clear. The radical 
conclusion (Brown 1995; Wettig & Brown 1996; Fryer & 
Kalogera 1997) is that the binaries in which neutron star 
enters a common envelope are not able to produce HMBP. 
This reduces the wealth of scenarios for the formation of 
HMBP to the following one: 

• 

II. (ms, ms) e , [gs, ms), (he, ms), {he, gs}, (he, he), (he, ns) e 
(he, ns) e , (ns, ns) e . 

In scenario II both binary components must initially have 
almost equal mass (within ~ 4%) and a primary with a 
relatively shallow convective envelope. This allows almost 
conservative mass exchange in the [gs, ms) phase. Due to 
the large amount of mass which is accumulated by the 
secondary star its evolution raps-up considerably and it 
can overfill its Roche lobe before the primary star explodes 
as a type lb supernova. The first-born ns is subsequently 
recycled in the stellar wind of its companion helium star. 
The stellar wind prevents the helium common envelope 
phase. As the outcome of hypercritical accretion is yet 
not clear, we consider below two options: the formation of 
a single black hole or the formation of a (bh, ns) binary 
for the case where both components do not merge within 
the common envelope formalism. 

Finally, there exist "wide" binaries which avoid mass 
exchange altogether. Both components may evolve inde- 
pendently, if the initial orbital period is <; 10 4 days or 
if the initial mass of the components exceeds 30-40 M Q 
and stellar wind drives the starts apart and prevents 
Roche-lobe overflow. The initial eccentricity is not af- 
fected by tidal effects and it is conserved until the first 
supernova event occurs. The eccentricity increases the sur- 
vival probability for the binary when the primary and the 
secondary subsequently explode in supernovae. It is likely 
that kicks prevent formation of (ns, ns) e pairs for (most of 
the) wide binaries. If one or two are formed anyway, both 
pulsars will be young and the binary will be visible as two 
single pulsars which happen to be spatially close for about 
10 Myr. Detection of binarity of pulsars is possible for or- 
bital periods below - 20 000 yr (Lamb & Lamb 1976). 



Binarity might be confirmed for very wide binaries if the 
proper motion of two stars are the same within about a 
km s^ 1 (Latham et al. 1984). Detection of a young pulsar 
in a wide binary may be an experimentum crucis for the 
hypothesis claiming that single stars and components of 
wide binaries do not produce radio pulsars because of the 
slow spin of the pre-supernova stars (Tutukov et al. 1984; 
Ibcn & Tutukov 1996). It will argue even against the oc- 
currence of kicks within a few tens of km s^ 1 (Portegies 
Zwart et al. 1997a). 

The neutron star in the close binary which was not dis- 
rupted by the first supernova cannot be detected as a radio 
pulsar due to free-free absorption of radio emission in the 
wind of its companion. Within 10 Myr to 20 Myr, before 
the second supernova occurs, the pulsar passes away. This 
first-born "dead" pulsar may be posthumously recycled 
in the stellar wind of its companion and/or during the 
common envelope phase (Bisnovatyi-Kogan & Romberg 
1975; Smarr & Blandford 1976; Bhattacharya & van den 
Hcuvel 1991). The interaction in a binary system is hy- 
pothesized to result in a sharp decrease of the magnetic 
field strength down to <~ 10 9 G (sec Bhattacharya 1996 for 
a list of suggested hypotheses and their classification). A 
weak field results in a prolonged lifetime for the recycled 
pulsar, about two orders of magnitude longer than for a 
"regular" pulsar born with a high magnetic field and with 
the same pulse period (see e.g. Fig. 1 in van den Heuvel 
& Lorimer 1996). 

2.2. Common envelopes 

The crucial stage in the formation of a HMBP in scenario I 
is the common envelope stage in which the envelope of the 
donor star is expelled from the binary system through the 
interaction with its companion. The net result of the com- 
mon envelope phase is a dramatic reduction in the orbital 
separation. Estimates of this reduction may be obtained 
from a comparison of the binding energy of the common 
envelope AE^ with the change of the orbital energy of the 
binary AE or b (see e.g. Iben & Livio 1993; Livio 1996 for 
discussion). The measure of this ratio is the common enve- 
lope parameter a cc = AE\, / AE OI b . Tutukov & Yungelson 
(1979) suggested the following equation for the variation 
of the separation of binary components in the common en- 
velope (the donor mass is M and its companion has mass 
to): 

(Mi + to) (M; — Mr) / 1 1 \ 
^ =«ceM f m(^---j, (1) 

where subscripts i and / refer to the initial and final states, 
respectively. This equation is based on the assumption 
that the spiral-in starts when the envelope surrounding 
the two cores has dimension of <~ 2a; and the envelope is 
ejected from the gravitational potential of both compo- 
nents. Recent 3-D smoothed-particle hydrodynamics cal- 
culations for the initial stage of formation of common en- 
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velopes (Rasio & Livio 1996) confirm such an assumption, 
at least qualitatively, by showing that very little mass is 
expelled during this stage and that the envelope expands 
considerably. 

On the other hand, Webbink (1984) suggested a more 
straightforward equation for otf ja\ variation: 



Mi (Mi - M f ) 



■q cc m 



M f 
2a f 



M 
2fli 



(2) 



where A is a structural constant which depends on the 
density profile in the star [A = 3/(5 — n) for polytropical 
models with index n; the usual assumption is A — 0.5], r/ co 
is the common envelope parameter and tl is the fractional 
radius of the Roche-lobe of the donor. Actually we handle 
the product ry ce A cXS (X single parameter. 

For completely non-conservative common envelopes, 
both equations result in equal Of/a; for a co w 4?7 co A. Con- 
frontation of population-synthesis calculations with ob- 
servations does not provide very rigorous constraints on 
the common envelope parameter. Nevertheless, popula- 
tion synthesis of galactic symbiotic stars (Yungelson et al. 
1995), Be/X-ray binaries (Portegies Zwart 1995), super- 
soft X-ray sources (Yungelson et al. 1996), close binary 
white dwarfs (Iben et al. 1997) and low-mass millisec- 
ond pulsars (Tauris 1996; Tauris & Bailes 1996) confirms 
that using a high value of the common-envelope param- 
eter gives a better agreement with the observations. Van 
den Hcuvel (1994) finds a reasonable evolutionary scenario 
with rj cc X = 2 for PSR J2145-0750, a low mass binary pul- 
sar, while for 77 ce A = 0.5 the spin-up of the neutron star in 
this particular case is problematic. Increasing ?7 C0 A to val- 
ues larger than unity decreases the discrepancy between 
the predicted and observed ratios of the birthrates for 
neutron star and black hole low-mass X-ray binaries (see 
Fig. 4 in Portegies Zwart et al. 1997b). On the basis of 3- 
D smoothed-particle hydrodynamics calculations Rasio & 
Livio (1996) claim that the deposition of energy into the 
common envelope must be highly efficient. Taken at face 
value, for Eq. 2 values of i] cc > 1 mean in that additional 
energy for the expulsion of the common envelope must 
come from sources other than the orbital energy. How- 
ever, Eqs. 1 and 2 are by no means anything more than 
indicative order of magnitude estimates. We argue below 
that reproduction of the distribution function for: the or- 
bital period, the eccentricity and the relative birthrate of 
observed HMBP are reconciable with high values for the 
common envelope parameter. 

2.3. The model 

For the present study we used the population synthe- 
sis code SeBa described by Portegies Zwart & Verbunt 
(1996). The basic difference with the Tutukov & Yun- 
gelson (1993a, b) code is that they use Eq. 1 instead of 
Eq. 2 which is used by Portegies Zwart & Verbunt for the 
treatment of binaries which evolve non-conservatively in 



the semi-detached and common-envelope phases. A de- 
tailed comparison of evolutionary sequences for models 
with similar initial parameters (Portegies Zwart & Ver- 
bunt 1996) shows that the sets of assumptions in the two 
codes result in reasonably similar configurations prior to 
the second supernova explosion in the system (although 
the stages in between can be quite different). Moreover, 
the enhancement of the amount of angular momentum 
loss during non-conservative mass transfer in the present 
model [rjj — 6 instead of 77,7 = 3, sec Eq. (38) in Porte- 
gies Zwart & Verbunt (1996)], makes the results more 
comparable 2 . Comparison of synthesis studies with the 
observed population of Be/X-ray binaries supports this 
choice for loss of angular momentum (Portegies Zwart 
1995). For the present study the mass-loss law in the stel- 
lar winds from the helium stars was accepted after Langer 
(1989): M = 5 x 10~ 8 M 2 - 5 M Q yr" 1 , which is about twice 
the value used by Portegies Zwart & Verbunt (1996). 

After the initial conditions for each binary are selected 
it is evolved until the binary merges, is disrupted or a pair 
of neutron stars is formed. For the latter, further decay 
of the orbit is computed from linearized general relativity 
(Peters 1964). At each moment in time, the distributions 
of the orbital elements are computed assuming a contin- 
uous formation of new systems and taking into account 
the decay of the orbital elements of the present (ns, ns) e 
binary population. 

The zero-age parameters of the simulated binary pop- 
ulation are selected from independent distribution func- 
tions. The mass of the primary M is chosen from a mass 
function of the form: 



dN oc M~ 2 - 5 dM. 



(3) 



The mass of the secondary star m is selected from a dis- 
tribution function for the mass ratio q a = m/M : 



*(«„) a (1 + qof- 



(4) 



We explored the dependence of the model upon <j) in 
the range from -4 to +5 and for a "flat" distribution 
ty(q ) = 1- For the semi-major axis distribution we ex- 
plored two options which are consistent with observational 
data: flat in logo between 0.1 R and 10 6 R (Abt 1983) 
and a Gaussian with a mean of 10 3 9 Rq and a dispersion 
of 1O 2 9 R0 (Duquennoy & Mayor 1991). The eccentricity 
for a binary was chosen between and 1 from the thermal 
distribution: S(e) = 2e. 

Observations of the space velocities of single radio pul- 
sars suggest that neutron stars may receive a kick at birth 
(Gunn & Ostriker 1970). A clear summing up of argu- 
ments for asymmetric kicks is given by van van den Heuvel 



2 For binaries with a mass ratio below ~ 0.4 assumption of 
rjj — 6 is similar to mass loss through the second Lagrangian 
point, for larger mass ratios less angular momentum is lost. 
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& Paradijs (1997) 3 . Following Hansen & Phinncy (1996) 
and Hartman (1997) we used the distribution for isotropic 
kick velocities (Paczyhski 1990) 



with u = v/cr and a = 600 km s -1 . Implementation of 
this distribution of kicks by Hansen & Phinney and Hart- 
man provides a reasonable model for the population of sin- 
gle pulsars close to Sun. As an alternatives we also made 
computations without kicks and with a Maxwellian kick- 
velocity distribution. 

All stars are supposed to be born in binaries. The to- 
tal birthrate is normalized to the present astration rate 
of 4 M Q yr~\ together with M min = 0.1 M and Eq. 3, 
which is consistent with observational estimates for the 
Galactic disc (see van den Hoek 1997; van den Hoek 
& de Jong 1997). The inaccuracy in the number of bi- 
naries formed with Mi > 8M Q due to uncertainties in 
the astration rate and the shape of mass-function below 
<~ 0.3M Q may be 0(2). The assumed age of the Galactic 
disc is T = 10 Gyr (Meynet et al. 1993; Carraro & Chiosi 
1994). The star formation rate is assumed to be constant 
throughout the evolution of the Galaxy. For all models 
a total of 5 x 10 5 binaries with a primary mass between 
8M Q and 100 M Q arc initialized. 

3. Results 

Our main purpose is to model the present Galactic pop- 
ulation of (ns, ns) e binaries with the aim to reproduce 
their birthrate, period-eccentricity distribution and to ob- 
tain an estimate for their merger rate and extrapolate it 
to the local Universe. 

3.1. Birth rates of neutron stars and related binaries 

In each of the eight models, for which we present the re- 
sults below, different assumptions are made concerning the 
kick velocity, the common envelope parameter, the mass 
ratio and semi-major axis distributions. Basic peculiarities 
of the models arc summarized in Table 1 . As a "standard" 
set of assumptions we consider the one with a flat initial 
distributions in q Q and log a, common envelope parameter 
r) cc \ = 2, and the kick distribution given by Eq. 5. 

In model A no kicks are imparted to neutron stars or 
black holes. 

Model B is identical to model A, but now neutron 
stars do receive a velocity kick upon birth chosen ran- 
domly from the distribution given by Eq. 5. This model 
is "standard" . Black holes also receive a kick, but scaled 
by mass: «bh = v ns m ns /mbh- No mass is ejected in the 
collapses to black holes, contrary to neutron stars. Small, 

3 See, however, Iben & Tutukov (1996) for arguments contra 
kicks. 



but not vanishingly, kicks for black holes may be consistent 
with the velocity dispersion of ~ 40 km s _1 of X-ray tran- 
sients with black-hole candidates (White & van Paradijs 
1996) which is higher than velocity dispersion of O-stars 
(~ 10 - 20km s" 1 , Wielen 1992). 

Model C is the only model in which the kick veloc- 
ity is selected from a Maxwellian distribution with a 3- 
dimensional velocity dispersion of a = 450 km s _ . 

Model D uses a smaller value for the common envelope 
parameter of 7y cc A = 0.5. 

The initial semi-major axes of orbits in model E are 
selected from a Gaussian (see sect. 2.3). 

In the models F and G the mass-ratio distribution is 
not flat: (f> = — 4 in model E and +5 in model G (see 
Eq. 4). 

Finally, in model H which used the standard assump- 
tions of model B, the accretion rate onto a ns in a common 
envelope is 10 8 x MEdd- All neutron stars which are able to 
increase their mass to > 2 M Q are to become black holes. 
No kicks are applied in these cases. 

Table 1 gives an overview of a selection of the com- 
puted models. The sum of the entries in the columns (2) 
to (5), and (8) gives for each model the annual rate of 
explosive events resulting in the birth of a young neutron 
star or a black hole. The estimated rate of formation of 
single neutron stars is given by the sum of the entries in 
columns (2), (6) and (7). In all models these rates are close 
to 1.5 x 10~ 2 yr _1 and are almost independent of the as- 
sumed semi-major axes and mass-ratio distributions of the 
binaries and the rate of binarity. This can be understood 
as a consequence of the relatively small number of high- 
mass secondaries. A smaller binary fraction together with 
the same astration rate results in an increase of the num- 
ber of single stars and therefore a higher birthrate for sin- 
gle pulsars. If all stars are born single the pulsar formation 
rate becomes 1.9 x 10~ 2 yr _1 . The model birthrate of pul- 
sars may be compared to the inferred Galactic occurrence 
rate of type II and Ib/c supernovae: (12 ± 6) x 10 _3 yr _1 
for SN II and (2 ± 1) x 10~ 3 yr" 1 for SN Ib/c (Cappellaro 
et al. 1997) 4 . Based on the modeling of population of pul- 
sars close to Sun, Hartman et al. (1997a) derive a pulsar 
formation rate of (9 ± 3) x 10~ 3 yr _1 , where a beaming 
factor of 0.3 was assumed. 

The relative birthrate of single and binary neutron 
stars is most strongly affected by the assumptions concern- 
ing kicks (compare model A with the models B to H, see 
Table 1). The formation rate of single non-recycled neu- 
tron stars, column (2) in Table 1, increases considerably 
when kicks are applied. This is not surprising because the 
majority of the wide non-interacting binaries and a consid- 
erable fraction of the binaries with shorter orbital period 
that experience a rather stable phase of mass transfer, 
which causes the orbital separation to increase, are disso- 

4 For the morphological type Sb, blue luminosity of 2 x 
10 10 L Q and H a = 75 km s" 1 Mpc" 1 . 
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Table 1. Results of the computations for a selection of the computed models. Columns (2) to (5) give the birthrates of single 
young radio pulsars, binaries which contain a young neutron star or a black hole and a non-compact star, (ns, *) and (bh, *), 
respectively, and binaries which contain a black hole and a young neutron star (bh, ns). Column (6) gives the birthrate of single 
recycled pulsars, column (7) - the occurrence rate of mergers between neutron stars and non-compact companions {ns, *} and 
column (8) gives the birthrate of HMBPs: (ns, ns). Columns (9) and (10) give the merger rate of neutron star binaries {ns, 
ns} and binaries that contain a neutron star and a black hole {bh, ns}, after 10 10 years of the evolution of the Galaxy with a 
constant star formation rate. Birthrates are normalized to galactic astration rate of 4 Mq yr -1 and 100% binarity. The numbers 
in parentheses for model H give the rates for the case when hypercritical accretion does not always result in the merger of the 
components. The statistical (Poissonian) uncertainty in the rates is smaller than 1%. 
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ciated upon the first supernova. On the other hand, the 
proportion of binaries which survive the first supernova 
explosion depends only weakly on the initial distributions 
for the mass ratio, semi-major axis, and the shape of the 
kick velocity distribution. The reason is that under any 
reasonable ^(q ) and 4>(a) most of the close binary sys- 
tems evolve to shorter periods and then it is mainly the 
average Ukick which determines the survival probability in 
the supernova explosions (see e.g. Fig. 8 in Lipunov et al. 
1996). 

The formation rate of bound neutron stars, (ns, *), 
column (3), drops by an order of magnitude when kicks 
are present. This entry in Table 1 gives a rough estimate 
of the formation rate of X-ray binaries. Since the majority 
of them are high-mass binaries, their total number in the 
Galaxy can be estimated. Actually, these systems gener- 
ally start as a Be/X-ray binary to evolve into an ordinary 
high- mass X-ray binary at a later instant. The X-ray life- 
time of a Be/X-ray binary is between 1 and 10 Myr (the 
remaining main-sequence lifetime of the B-star compan- 
ion of the neutron star). The estimated galactic popu- 
lation is consequently somewhere between 10 3 and 10 4 . 
The lifetime of high-mass X-ray binaries is of the order 
of 10 3 yr to 10 5 yr, depending on the evolutionary status 
of the donor (e.g. Savonije 1979; Massevich et al. 1979; 
Hcllings & DeLoore 1986). This results in a total of sev- 
eral tens of high-mass X-ray binaries in the Galaxy. These 
estimates of numbers of Be- and high-mass X-ray sources 
are consistent with observational estimates of > 2000 and 
55 ± 27, respectively (Meurs & van den Heuvel 1989, and 
50 to 80 by Iben et al. 1995). Models without a kick seem 



to overestimate the numbers of X-ray binaries by an order 
of magnitude. 

Binaries which contain a black hole and a stellar com- 
panion are not as common as binaries which contain a 
neutron star. If their lifetime is comparable to the life- 
time of high-mass X-ray binaries which contain a neutron 
star, their number has to be ~ 1/10 of the former, in a 
agreement with the observations. 

Derivation of the expected birthrate and the number 
of low-mass X-ray binaries from our calculations would 
require a more rigorous consideration of very severe con- 
straints on the orbital periods and masses of donors in 
prospective candidates (Kalogera & Webbink 1996) and 
is beyond the scope of the present paper. 

Entries in columns (6) to (8) of Table 1 refer to neu- 
tron stars which were born first in the binaries and were 
subsequently recycled or merged with their stellar com- 
panions. Evidently, the interaction with companions in- 
fluences the radio properties of neutron stars. We still call 
them "recycled" pulsars, though decay of the magnetic 
field and/or the action of the propeller mechanism can 
prevent a certain fraction of them to show up in radio. The 
birthrate of single recycled pulsars is given in column (6) 
of Table 1. Observational properties of these objects are 
discussed in detail by Hartman et al. (1997b). Typically, 
single recycled pulsars contribute for less than 1% to the 
total population of single radio pulsars 5 . This is consistent 
with the presence of eight isolated objects among observed 

5 The birthrates found by Hartman et al. are smaller by 
about a factor two because they assumed a binarity fraction of 
50%. 
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millisecond pulsars in the Galactic disc. In the absence of 
kicks (model A) <~ 1/7 of all pulsars are old and possibly 
recycled; this contradicts the observations (see however 
Deshpande et al. 1995 for arguments which suggest a high 
fraction of recycled pulsars) . 

The presence of kicks significantly increases the pro- 
portion of systems in which the neutron star merges with 
its stellar companion. The main reason for this is that 
a kick can decrease the orbital period or even shoot the 
neutron star directly into its companion (see also Brandt 
& Podsiadlowski 1995). In the first case, the smaller or- 
bital separation of the binary causes a merger to occur 
more easily after a phase of unstable mass transfer. If the 
merger does not lead to the formation of a black hole the 
neutron star may, after ejection of the common envelope, 
emerge as a single recycled pulsar. These systems may in- 
crease the fraction of recycled pulsars by a factor of two. 

3.2. Neutron star-neutron star binaries 

Comparison of column (3) with column (8) from Table 1 
shows that less than ~ 5% of the Be- or high-mass X- 
ray binaries survive the second supernova and produce a 
recycled binary pulsar, about 25% is dissociated in the 
second supernova and about 10% merge into a single ob- 
ject. The remaining ~ 60% never experience the second 
supernova. A minor proportion of the latter systems may 
become LMXB and, later, millisecond pulsars. 

We confirm the result from previous computations that 
the birthrate of (ns, ns) binaries drops by more than an 
order of magnitude when kicks are implied. This reflects 
in the birthrates in Table 1 . More binaries are dissociated 
when the kick velocity is taken from a Maxwcllian distri- 
bution instead of the distribution of Eq. 5, because of the 
larger average magnitude of the kick. 

Directly after the second supernova, a binary contain- 
ing two neutron stars of which one is visible as a young ra- 
dio pulsar and the other may be recycled. The later born 
pulsar dies after about 10 Myr. A recycled pulsar may 
live up to a few billion years. For most of this time span 
the binary may be detectable as a recycled pulsar with 
a "dead" companion. A few known systems are listed in 
Table 2. For PSR J1518+49 and B2303+46 only the total 
mass is known and they may have a massive white dwarf 
instead of a neutron star as a companion. White dwarf 
companions will be ruled out if optical observations will 
not discover objects of a visual magnitude of <~ 26 for 
PSR J1518+49 and between 23 to 25 for the other pulsar 
(Fryer & Kalogera 1997). 

The long pulse period and strong magnetic field of 
PSR B2303+46 do not exclude that it is the second born 
pulsar which we observe. In this case, its companion would 
be at most ~2x 10 7 yr older (the evolutionary timescale 
of a M <~ 10 M Q star). The companion is cither a dead 
pulsar which did not resurrect as a recycled pulsar or it 
is beaming away from our planet. This discussion also ap- 



plies to PSR B1820-11 for which even the total mass is not 
known. It might be the young pulsar with a main-sequence 
star, a white dwarf or an unseen pulsar as companion (see 
Phinney & Verbunt 1991 and Portegies Zwart & Yungel- 
son (1997), for discussion). 

The birthrate of radio pulsars in neutron star bina- 
ries must equal the birthrate of recycled pulsars in similar 
systems (Bailes 1996). The absence of "young" pulsars 
in pairs with old ones among ~ 650 objects restricts the 
relative birthrate of (ns, ns) systems to 1/650. Bailes' cri- 
terion does not restrict the absolute birthrate of (ns, ns) 
binaries; it lacks absolute calibration. Except for model A, 
none of the models contradicts this statement. If it is ac- 
tually the young pulsar we observe in PSR B2303+46 and 
B1820-11 the limit on the relative birthrate increases to 
2/650. 

3.3. Orbital period-eccentricity plane for neutron-star bi- 
naries 

The only observational parameters of HMBPs to which 
one may compare the results of evolutionary computations 
and population synthesis are the orbital period P or b and 
eccentricity e. Other characteristics of the binary like the 
center-of-mass velocity are ill known from observations. 
For a comparison of the characteristics of the individual 
pulsar, like the strength of the magnetic field, the life- 
time and the pulse period, no theory is available. In the 
Porb ~ e diagram it is only the emission of gravitational 
waves that effects the further evolution of the orbital pa- 
rameters. 

Figure 1 gives the distribution function for the orbital 
parameters of (ns, ns) binaries at birth for a model with- 
out a kick (A) and for model B with a kick velocity dis- 
tribution according to Eq. 5. The population of (ns, ns) 
binaries with periods beyond 10 4 days is not shown. Since 
in models without kick most wide binaries survive the su- 
pernovae explosions bound, ~ 70% of all (ns,ns) have or- 
bital periods greater than 10 4 day. Nondetection of young 
pulsars in wide binaries then argues for kicks of at least of 
several tens of km s~ , see Portegies Zwart et al. 1997a. 
In the models with a kick this population contributes only 
little 10%, see also Portegies Zwart & Verbunt 1996). 
In both models the adopted common envelope parameter 
in Eq. 2 is 77 CC A = 2 (which is roughly equivalent to a cc = 
0.5 in Eq. 1). The reason for this choice is the inability 
for smaller ?7 CC A to reproduce the position of the observed 
high-mass recycled pulsars in the P or b - e plane. Partic- 
ularly, the probability for the formation of binary pulsars 
with P or b £ 2 days and e ^ 0.3 and P or b ~ 10 days and 
e <~ 0.6 — 0.7 appears to be very small. The reason is that 
low efficiency of the deposition of orbital energy into the 
common envelope results in the formation of very close bi- 
naries. For comparison we give the P or b - e diagram for 
??ccA = 0.5 in Fig. 4. 
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Table 2. Observed population of massive neutron star binaries. Their names are followed by the pulse period, the period 
derivative and the orbital elements, age and magnetic field strength. References: 1 - Nice et al. (1996), 2 - Wolszczan (1990), 3 
- Lyne & McKenna (1989), Bhattacharya & van den Heuvel 1991, 4 - Hulse & Taylor (1975), 5 - Deich & Kulkarni (1996), 6 - 
Stokes et al. (1985). Data which is not available in the cited papers were taken from Taylor et al. (1993). 
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Fig. 1. The probability distribution for the orbital parame- 
ters of the (ns, ns) binaries with a period smaller than 10 4 
days at the moment of birth. The upper panel gives the prob- 
ability distribution for the model with n cc \ = 2, flat initial 
distributions over q a and log a and without a kick; lower panel 
is for the same model with a distribution of kicks according 
to Eq. 5. Darker shades indicate higher probability to find a 
neutron star binary at that location in the period-eccentricity 
plane. The gray scaling is given in the upper left corner of each 
figure, the darkest shade in the upper panel corresponds to a 
birthrate of 1.2 x 10~ 5 yr~\ for the lower frame this value is 
3.5 x 10~ 7 yr _1 . The integrated birthrates are given in Table 1, 
column (8) 



If no kicks are imparted the P or b ~ e diagram has 
clear-cut borders (see Fig. 1). The shortest orbital pe- 
riod is determined by the requirement to accommodate 
a helium star with a neutron star companion in the pre- 
explosion orbit. The supernova tends to widen the binary 
orbit and the minimum eccentricity is obtained from the 
fact that at the instant of explosion all pre-neutron stars 
still have helium envelopes of AM >, 0.5M Q . A minor pop- 
ulation with low e and 20 £ P or b ~ 30 day descends from 
massive systems which had eccentric orbits before the 
second explosion. The major, short-period, population in 
Fig. 1 descends from systems which evolved via scenario I 
with all mass exchanges being unstable, while wider sys- 
tems had at least one phase of stable mass exchange be- 
fore the formation of the first neutron star. Presence of 
these two populations was also clearly seen in the figures 
in Tutukov & Yungelson (1993a) and Portegies Zwart & 
Spreeuw (1996) papers. The lower panel in Fig. 1 demon- 
strates that, apart from reducing the birthrate, kicks have 
the effect of "smearing" the period-eccentricity distribu- 
tion of the neutron star pairs. The eccentricity distribution 
is more strongly affected by the presence of a kick than the 
period distribution. Striking in the lower panel in Fig. 1 is 
the existence of a population of binary neutron stars with 
orbital periods below ~ 0.1 day, in contrast to the compu- 
tations without a kick. These short periods are the result 
of a kick in the "right" direction to reduce the orbital 
period upon the second supernova. Due to the emission 
of gravitational waves this population is extremely short- 
living. If kicks are absent, only the Hulse- Taylor pulsar is 
born in the region in the P or t, - e diagram with high for- 
mation probability. The extremely small probability for 
the formation of zero-age pulsars with parameters similar 
to that of J1518+49 and B2303+46 is appealing. 

Figure 2 shows the effect of time evolution on the pop- 
ulation of neutron star binaries from the standard model B 
with a kick, giving the probability distribution of them at 
the age of 100 Myr, 1 Gyr and 10 Gyr. The limit of 100 Myr 
is suggested by the characteristic age of PSR B1913+16. 
Figure 2 may be compared with the lower panel of Fig. 1 



Simon Portegies Zwart & Lev Yungelson: Formation and Evolution of Neutron Star Binaries 



9 



which gives the probability distribution of these binaries at 
birth (which is more appropriate for relatively young pul- 
sars like PSR B2303+46). PSR B1913+16 fits extremely 
well into the model population for its age. The rest of the 
object fall into regions with probability by a factor of a 
few lower, but still relatively high. The observed system 
may be found in the short period area of the diagram 
if the birthrate of its progenitors is high or if the pulsar 
is old. The first is probably the case for PSR B1913+16 
and the other option is for older PSR B1534+12. Pulsar 
PSR J1518+49 is found in a low-birthrate area of P or b - e 
diagram which suggests a long age. Relatively young pul- 
sars, like PSR B2303+46, may be found in any area of the 
Porb - e diagram except for the very short periods where 
the merger time is smaller than the characteristic age. 

The probability distribution shifts to larger orbital 
period as the population becomes older, because short- 
period objects merge due to the influence of gravitational 
wave radiation while the long-period binaries survive for 
a longer time. As the age of the population increases pul- 
sars older than PRS B1913+16 are found in the area of 
the P orh - e diagram with a higher probability; the for- 
mer pulsar drifts to the lower probability area. The mid- 
dle panel of Fig. 2 fits best with the observed four systems 
suggesting that the average age of the observed popula- 
tion is close to 1 Gyr. It is, however, hard to quantify this 
statement. Particularly, there may be a differential selec- 
tion effects related to the different lifetimes of recycled 
pulsars, which depends on their specific accretion history. 

The bottom panel of Fig. 2 represents the total pop- 
ulation of binary neutron stars currently present in the 
Galaxy. The (ns, ns) binaries we observe today are a sub- 
set of this population. Since typical merger time for most 
(ns, ns) binaries is several hundred Myr (see Fig. 7 lower 
panel) and radio lifetime of recycled pulsars is ^ 1 Myr 
(van den Heuvel & Lorimer 1996), the majority of the 
neutron star binaries present in the Galaxy are dead! 

Figure 3 gives the probability distribution for (ns, ns) 
pairs younger than 100 Myr in the P or b - e plane for the 
model where the kicks are selected from a Maxwellian ve- 
locity distribution. The gray shading in this Figure is nor- 
malized to that in the upper panel of Fig. 2. This model 
has a smaller (ns, ns) birthrate than model B, this re- 
flects in the absence of the darkest shades in Fig. 3. The 
presence of a non-negligible, yet not observed, popula- 
tion of (ns, ns) binaries with orbital periods larger than 
~ 10 days in model B (see Fig. 2) reflects the higher prob- 
ability of low-velocity kicks in the distribution given by 
Eq. 5 relative to the Maxwellian. Though the latter model 
suggests almost equal probability for the detection for 
PSR B1534+12, B1913+16 and B2303+46, it gives an ex- 
tremely low birthrate for binaries with orbital parameters 
similar to PSR J1518+49. Even high age does not increase 
the model galactic population of pulsars with similar or- 
bital parameters. This discrepancy does not go away when 
the velocity dispersion of the Maxwellian distribution is 



decreased by a factor of 2. (Hansen & Phinney 1997, sug- 
gested that a Maxwellian kick-velocity distribution with 
a = 190 km s _1 fits satisfactory with the observed single 
pulsar velocities.) 

Model E, with initial Gaussian distribution over logo 
differs only slightly from the model with an initial flat 
distribution in log a. It gives a similar birthrate of (ns, 
ns) binaries and reproduces the P or b - e relation equally 
satisfactory. This means that most of (ns, ns) stars de- 
scend from initial systems having orbital separations in 
the interval where the probabilities for formation for both 
distributions are comparable. 

Our computations included several additional runs 
with different parameters which produced models which 
we consider less satisfactory than the standard model. 
Model F with initial ^-distribution strongly raising to 
((f) = —4, in Eq. 4) gives birthrates of single neutron 
stars and (ns, ns) binaries similar to these in the standard 
model B, but does not reproduce the P or b - e distribution 
satisfactory. Particularly, it does not reproduce the pop- 
ulation with P or b ~ 10 day and e £ 0.7 and pulsars with 
Porb ~ 0.4 days and e s=s 0.2 — 0.3. This is an understand- 
able result of the high proportion of mergers in the first 
(hydrogen) common envelope event and of strong reduc- 
tion of the orbital separation in the systems which avoid 
mergers. Model G, with q Q strongly raising to 1 (<j> = +5), 
does not basically differ from the other models in the total 
production rate of pulsars and it reproduces the P or b - e 
relation rather well. However, it is the least satisfactory 
model respective to the ratio of formation rates of binary 
and single pulsars (see discussion section). 

Figure 4 gives the P or b - e plane for the population 
younger than 100 Myr for (ns, ns) systems of model D. In 
this model the efficiency for the deposition of orbital en- 
ergy into the common envelope during the spiral-in is con- 
siderably smaller than in the standard model: fy cc A = 0.5. 
The birthrate of (ns, ns) binaries in this model is consid- 
erably smaller than in the standard model, as a result of 
greater proportion of mergers under low r/ cc . However, all 
four HMBPs appear in low-probability areas of the dia- 
gram. This was the basic reason for pushing rj cc to higher 
values in our efforts to obtain satisfactory models for the 
population of binary neutron stars. 

Figure 5 gives the P or b - e distribution of the neutron 
star binaries which originate from model H for two age lim- 
its. This distribution is hardly affected by age (compare 
with model B in Fig. 2). Since systems which experience 
a common envelope are absent in this model, the average 
orbital separation is relatively large. The distribution is 
static for systems with P or b <; 1 day; it is hard to discrim- 
inate between the 100 Myr and the 10 Gyr populations. 

The population which originates from model H is a 
subset of model B, and Fig. 5 consequently represents the 
P or b - e distribution of this subset in Fig. 2. Only the 
young pulsar PSR B2303+46 is present in the area with 
the highest probability in Fig. 5. Therefore it is appeal- 
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ing to suggest that it originates from scenario II . Simi- 
larly, van den Heuvel et al (1994) argue that the evolution 
of the progenitor of this system was dominated by mass 
loss in a stellar wind rather than by Roche-lobe overflow. 
A common envelope phase is required to form short pe- 
riod binary pulsars. The apparent "success" of Wettig & 
Brown (1996) and Fryer & Kalogera (1997) in reproduc- 
ing PSR B1913+16 and B1534+12 is mainly due to their 
fine tuning of the parameters of the pre-supernova binary. 
Such systems are not the dominant contributors to short 
orbital period HMBP. 

Figure 6 gives the cumulative velocity distribution for 
the center of mass of the merging (ns, ns) binaries which 
results from the model computations. For comparison, we 
also give the kick distribution (Eq. 5). As expected, the 
model without the kick produces low velocity binary neu- 
tron stars and models B and C produce the higher velocity 
ones. The low-velocity shoulder for the model without a 
kick (left dash-dotted line), shown for comparison, is the 
result of very wide - non-interacting - binaries which hap- 
pen to stay bound due to the orbital eccentricity and low 
mass loss during the supernova. The short period (ns, ns) 
binaries from the model without a kick reach much higher 
velocities (right dash-dotted line), but by far not as high 
as in the models where a kick is imparted to the newly 
formed neutron star. The difference between the models 
with a Maxwellian kick and a Paczyhski kick becomes 
quite noticeable in the tail of the distribution: model B 
produces more high velocity neutron star binaries. The 
velocity distribution for HMBP which are produced in 
model H is very similar to that of model B (solid line) 
but it is shifted to lower velocities by about 50km/s, as 
a result of the absence of very close systems among their 
progenitors. The velocity distribution for the single pul- 
sars is similar to the kick distribution, the contamination 
by recycled pulsars (which typically have smaller veloci- 
ties, Hartman et al. 1997b) is negligible. The distribution 
for HMBP has a lower proportion of high- velocity objects 
compared to the single pulsars because the high velocity 
tail of the kick distribution prevents the formation of (ns, 
ns) binaries. There are not references to the spatial veloc- 
ities of binary pulsars in the refereed literature and it is 
impossible to confront Fig. 6 with observations. We should 
note, however, that discovery of a HMBP with velocity in 
excess of <~ 250 km s _1 would argue for the presence of 
kicks. 

4. Discussion 

4-1. Populaton of HMBP and their merger times 

Population synthesis computations for the formation and 
further evolution of binary neutron stars provide estimates 
for their birthrate and the expected orbital elements of 
the present day population. The most crucial parameter 
for modeling of (ns, ns) binaries is the common envelope 



parameter. This parameter is not strongly constrained by 
observations. Population synthesis and smoothed parti- 
cles hydrodynamic computations suggest that energy de- 
posited into common envelopes may be high compared to 
the orbital binding energy. Our preferred model B is char- 
acterized by such a high efficiency: n cc A = 2. Higher values 
of ?7 C cA up to 4 still result in agreement with the observed 
population but the birthrate of neutron star binaries be- 
comes uncomfortably high compared to the formation rate 
of single radio pulsars. For n cc A ^ 1 the birthrate of neu- 
tron star binaries decreases significantly but the period- 
eccentricity diagram fits badly with the four observed sys- 
tems. However, it might well be that n co A depends on the 
evolutionary stage of the binary; for an inspiralling neu- 
tron star this parameter might have a very different value 
than for a main-sequence star or white dwarf. Neutron 
stars can be very efficient in dispelling the common en- 
velope by energy released by e.g. small explosions (Fryer 
et al. 1996) or semi-steady burning on its surface during 
the spiral- in process. 

Only the models in which a kick is imparted to the 
newly born neutron star from the distribution proposed 
by Paczyhski (1990) with parameters suggested by Hart- 
man et al. (1997a) give a satisfactory birthrate as well 
as a distribution over P or i, — e. Increasing the fraction of 
low-velocity kicks in Paczyhski distribution will however, 
increase the quality of the fit of the P or b - e plane. This 
distribution seems at the moment preferable because the 
reproduction of transverse velocities of young pulsars with 
measured proper motions requires the presence of a sig- 
nificant number of objects with low velocities (Hansen & 
Phinney 1996, Hartman et al. 1997a). The actual kick ve- 
locity distribution might have an even larger contribution 
of low velocity kicks than Eq. 5 suggests. 

Apart from kick velocity, all other parameters have 
a relatively small effect. Birth- and merger rates vary 
within a factor of a few, but not by orders of magnitude. 
Therefore it is hard to decide on population estimates of 
birthrates which model parameters are favoured. 

The observations may be best fit with the computed 
probability distribution in the P or b ~ e plane if the age 
of the recycled neutron stars is of the order of 100 mil- 
lion to 1 billion year. Note that there might well be ef- 
fects other than the age alone which have its influence 
on the detectability of neutron star binaries; the pulse 
period, the strength of the magnetic field and the lumi- 
nosity. All these parameters give rise to a different life- 
time of the pulsar and different observational character- 
istics which make the selection effects hard to quantify. 
It is appealing that the youngest system PSR B2303+46 
fits the relatively high probability region of P or b - e di- 
agram for the youngest pulsars (Fig. 1, lower panel), the 
oldest one PSR J1518+49 fits best if the age of the pop- 
ulation is as old as a few billion years (close to the one in 
the lower panel of Fig. 2) and intermediate age pulsars 
PSR B1913+16 and PSR B1534+12 are best reproduced 
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by a population of a few hundred Myr (Fig. 2). Since the 
best fit between the observations and the models is at an 
age of the population between 100 Myr and 1 Gyr, this 
suggests that the majority of the recycled neutron stars 
are born with low magnetic field. 

Taking into consideration the characteristic age, or- 
bital period and the eccentricity of PSR B2127+11C, one 
may suggest that it is formed by scenario I instead of 
an exchange interaction in the high-density core of the 
globular cluster M15 (see Phinney & Sigurdsson 1991). 
However, to be primordial, the pulsar should be at least 
10 Gyr old (the age of the globular cluster). This pulsar 
has to merge with its companion within a few 100 Myr 
and it is highly unlikely to observe it in the last percent 
of its lifetime as a (ns, ns) binary. 

The birthrate of binary pulsars relative to singles is 
~ 1/430 in the standard model, while in the most extreme 
model H it is <~ 1/1775. These numbers have to be reduced 
for the binary fraction, which is smaller than 100%. Ac- 
cording to van den Heuvel & Lorimer (1996) the lifetime 
of recycled pulsars is a factor ~ 3 larger than their char- 
acteristic age and therefore it is between ~ 100 Myr and 
~ 1 Gyr. On the other hand population synthesis studies 
suggest that ~ 90% of the recycled pulsars are missed in 
pulsar surveys (Curran & Lorimer 1995). The lifetime of 
young pulsars is ~ 10 Myr. Combining these numbers we 
arrive at the conclusion that the relative birthrate of sin- 
gle pulsars and binary pulsars also reflects their relative 
galactic number density. If beaming factors are similar for 
old and young pulsars, relative numbers of ns and (ns, ns) 
in our models do not contradict observations (see Bailes 
1996). As the birthrate of single pulsars in our model is 
consistent with the rate of supernovae and the birthrate 
of single pulsars, we may conclude that our models give 
basically the actual Galactic birthrate of binary neutron 
stars. 

A fraction of the neutron stars in binaries may 
never become observable as recycled pulsars. The large 
pulse periods of X-ray pulsars suggest that observed re- 
cycled pulsars are resurrected at later stages. It may 
be shown that in the common-envelope stage that fol- 
lows the X-ray pulsar stage, if B Z, 10 11 G, accre- 
tion at the Eddington rate doesn't spin neutron stars 
to P ^ 100 ms (Ergma & Yungelson 1977). Recycling 
of a pulsar to P cq <~ 30 ms requires the accretion of 
AM w 0.17(P cq /ms)~ 4/3 w 0.002M o (e.g. Bhattacharya 
1996). The only evolutionary stage in which this can hap- 
pen is during accretion from the wind of the accompany- 
ing helium star. The "propeller mechanism" (Illarionov & 
Sunyaev 1975), which prevents accretion onto the neutron 
star may cause the neutron star to remain dead. Ergma & 
Yungelson (1997) derive an estimate for the orbital period 
at which the propeller mechanism is effective: 

Porb £ 0.5M n 3 s /V 3 00 (M ns + M he )" 1 /2M h 1 e 5/8 hr, (6) 



where wiooo is the velocity of the helium-star wind in 
1000km s _1 and masses are in Mq. The spin slow-down 
time in the (ns, he) systems is comparable to the helium 
burning time. For typical wiooo ~ 1 - 2, the critical period 
is of the order of only several hours for a low mass helium 
star, which is quite short. If the propeller mechanism is ef- 
fective, accretion is possible for short period binaries with 
a low mass helium star or for wide systems with a high 
mass helium star. In the computations of Ergma & Yun- 
gelson it appeared that the majority of the neutron stars 
are not able to accrete from helium companion winds 6 . 
The first born pulsar has plenty of time (10 to 20 Myr) 
to die before its companion becomes a neutron star. If 
the propeller mechanism is effective, the first born pul- 
sar will never be resurrected from the graveyard. In this 
case the birthrate of recycled pulsars in binaries may be 
smaller than the birthrate of neutron star binaries. Ac- 
cording to Lipunov et al. (1997) this fraction may be as 
small as ~ 1%, depending on the assumed initial mag- 
netic moment of the neutron star and the decay timescale 
of its magnetic field. There may be one such a binary 
among the observed pulsars, and in fact it is possible that 
PSR B2303+46 contains a non-revived pulsar next to a 
young one (see Table 2). 

An observational estimate for the birthrate of HMBP 
can be derived from the estimate of their number in the 
Galaxy and their lifetimes as pulsars (Kulkarni & Narayan 
1988). Using scaling factors, correction factors for incom- 
pleteness of the sample and beaming factor from Cur- 
ran & Lorimer (1995) and lifetimes from van den Heuvel 
& Lorimer (1996) for PSR B1534+12, B1913+16, and 
B2303+46 we arrive at a birthrate of 2.7 x 10 -5 yr -1 . 
This is comfortably close to our standard model B (Ta- 
ble 1). This derived birthrate heavily relies on the young 
PSRB2303+46. Excluding it reduces the birthrate to 
0.8 x 10 _5 yr _1 , which is obviously the estimate obtained 
by van den Heuvel & Lorimer (1996). There is no data 
available on scaling factor for PSR J1518+49, but we may 
expect that it will enter the estimate of the birthrate with 
a very small weight because of its long characteristic age. 

To summarize, we may argue that our standard model 
B, which reproduces best the orbital periods and eccen- 
tricities of the observed HMBP is not in conflict with ob- 
servations also in respect to the birthrate of these systems. 
Therefore we may expect, that this model gives a realis- 
tic estimate of the merger rate of neutron star binaries 
w 2 x 10~ 5 yr _1 . This number is not inconsistent with the 
estimated rate of 0.8 x 10~ 5 yr _1 from van den Heuvel & 
Lorimer (1996) based on two observed pulsars with merger 
times shorter than the Hubble time and the assumption 

6 The results of Ergma & Yungelson cannot be directly ap- 
plied here for numerical estimates, since their population syn- 
thesis was made under different assumptions. Qualitatively it 
is clear that for the present computations accretion would be 
prohibited in even larger proportion of all systems, since here 
we typically produce wider (ns, he) binaries. 
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of a steady state in which the birthrate equals the merger 
rate. Actually, from Table 1 and Figs. 1 and 7 it is imme- 
diately clear that there does not exist a steady state in 
which the birthrate equals the merger rate. About 10% to 
20% of all (ns, ns) binaries merge after the recycled pul- 
sar has died within 0.5 - 1 Gyr (see Fig 7). These binaries 
hang around undetected in the Galaxy until they merge. A 
considerable fraction of these (20% to 30%) are parked in 
wide orbits in the Galactic halo or escape from its poten- 
tial well (sec Fig. 8 and 9). Another ~ 10% merge within 
~ 10 Myr, before even the young pulsar dies (see Fig 7). 
These systems may also avoid detection due to their short 
lifetime (their number relative to the number of single 
pulsars is only ~ 1/4000 in model B for 100% binarity). 
Figure 7 suggests a revision of the "observational" esti- 
mate upward by a factor -~l.5tol.2xl0 -5 . This correc- 
tion factor my increase even further taking non-recycled, 
hence non-observable, close binary neutron stars into ac- 
count. (Curiously enough the situation with merging neu- 
tron stars is similar to that of merging white dwarfs which 
are candidate progenitors for SNe la; most of the mergers 
will occur when both white dwarfs have cooled beyond 
detectability, Iben et al. 1997.) 

On the other hand, about 1/3 of all (ns, ns) binaries 
born annually in the Galaxy have merger times exceeding 
10 Gyr and are accumulating in it. The total number of 
dead non-merging pairs present in the Galactic disc may 
be as large as 10 5 . 

Figure 7 indicates that the Galactic disc merger rate 
is determined mainly by the star formation rate during 
the last 1-2 Gyr. The average star formation rate over the 
history of the Milky Way is about 4 times higher than the 
current rate (e.g. van den Hoek 1997). The same is the case 
for other massive (~ 10 11 M ) spiral and irregular galax- 
ies (Gallagher et al. 1984; Sandage 1986). The star forma- 
tion rate was much higher than present only in the first few 
Gyr of the existence of the Galaxy. Figure 7 shows that 
this early star burst contributes only little to the current 
merger rate of (ns, ns) binaries. The same point has to be 
made with respect to elliptical galaxies. As most stars in 
them were formed in a short ( ^ 1 Gyr) burst, overwhelm- 
ing majority of merger candidates already merged long 
ago. Figure 7 suggests that in a galaxy which was formed 
in a 1 Gyr burst, the present rate of (ns, ns) mergers per 
unit mass is ^ (10 — 15)% of that in the galaxy of the 
same mass but formed in 10 Gyr. In the extrapolation of 
the Galactic rate of merger to the local Universe the dif- 
ferent contributions from spirals and ellipticals have to be 
considered (see Bagot 1996 for a parametric study of this 
issue). Scaling to the total number of galaxies (spirals as 
well as ellipticals) within 100 Mpc after Curran & Lorimer 
(1995) gives an upper limit for the number of mergers 
detectable by the advanced LIGO detector of lyr -1 for 
model B. 



4-2. Merger rate in the Universe andj-ray bursts 

A rough estimate of the total number of merger events in 
the Universe may be obtained from scaling the Galactic 
merger rate by the ratio of blue-band luminosities of the 
Milky Way and the Universe: R s=s 10 -2 /iMpc -3 , where 
h is Hubble constant in 100 km s _1 (Phinney 1991). For 
H = 75 km s~ 1 Mpc" 1 this gives, extrapolating model B, 
~ 4.6 x 10 events per year. This number may be compared 
to the number of 7-ray bursts, which are also suggested 
to originate from neutron star mergers (Paczyhski 1986). 
The observational estimate of the rate of detectable 7-ray 
bursts in the Universe is several times 10 -6 per year per 
galaxy (Mao & Paczyhski 1992; Piran 1996). The discrep- 
ancy with the observed rate (for model B) may be erased 
if 7-ray emission is beamed into an opening angle of ~ 10° 
(Mao & Yi 1994; Piran 1996). The anisotropy of the en- 
ergy release in 7-ray production is suggested by numerical 
merger models (Davies et al. 1994). 

Virtually all (ns, ns) binaries in the models with a 
kick escape from the potential well of a globular star clus- 
ter (see Fig. 6) and a considerable fraction is expected 
to escape from a dwarf galaxy. In the model without a 
kick, this fraction is only ~ 25% (about 10% for dwarf 
galaxies). But in the latter model the majority of (ns, 
ns) systems originate from wide non-interacting binaries 
which never produce a recycled pulsar and which will not 
merge within the age of the Universe. The (ns, ns) pairs 
confined to parental stellar systems have typically large 
orbital periods and are therefore very fragile for three- 
body encounters. The result is that globular clusters are 
deficient of HMBP's according to the models with a kick 
and also according to the models without a kick! In the 
presence of a kick the formation of (ns, ns) binaries is not 
an efficient way to keep a neutron star bound to a globular 
cluster. 

The maximum distance traveled by a (ns, ns) binary 
before it merges may be computed using the distribution 
functions for the merger time and the center-of-mass ve- 
locity of the binaries. Figure 8 provides this distribution 
for the HMBPs. It is interesting to notice that in the mod- 
els without a kick the typical distance traveled is consid- 
erably larger than in the models with a kick. Note also 
that the distributions for both models with a kick are in- 
distinguishable. The binaries in the high-velocity tail of 
model B (see Fig. 6) merge on a short time scale (see 
Fig. 7) whereas the low- velocity binaries in the model with 
a Maxwcllian kick live rather long. 

The distance traveled perpendicular to the galactic 
plane in the potential of the Galaxy can be estimated from 
a relation between the z-componcnt of the initial velocity 
v(z) of an object and the maximum galactic height z (Tu- 
tukov & Yungelson 1978): 

v(z) w [6z - 2000 ln(l + 0.003z)] 1/2 , (7) 
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where v is in km s~ and z in pc. The maximum z 
reached by a (ns, ns) binary before coalescence was es- 
timated using the velocity from the distribution given in 
Fig. 6 divided by \f?> (assuming that the velocity is dis- 
tributed isotropically) . To correct for the pre-merger life- 
times which may be shorter than the travel time to the 
maximum z, the minimum of z resulting from Eq. 7 and 
the maximum distance traveled (Fig. 8) was adopted. 

Figure. 9 shows that for the models without a kick 
the majority of the systems do not reach a high z (left 
dash-dotted line), as expected from the enormous low- 
velocity tail in Fig. 6. The short period (ns, ns) binaries 
from model A reach considerably higher Galactic heights 
(right dash-dotted line) but by far not as high as in the 
models with a kick. Note that the majority of the binaries 
in the low-z shoulder (left dash-dotted line in Fig. 9) are 
not visible as systems with a recycled pulsar, since they 
did not experience a phase of mass transfer. The (ns, ns) 
binaries from the models where a kick is imparted to the 
newly formed neutron star show a much more extended z 
distribution. Actually, a considerable fraction of these sys- 
tems escape from the potential of the Galaxy or are parked 
in the Galactic halo. Escape of some (ns, ns) beyond the 
sensitivity limits of detectors may be a factor contributing 
to the underestimate of the galactic birthrate of HMBP. 
Since the life times of single radio pulsars are considerably 
shorter than that of recycled pulsars, of which the lifetime 
is comparable to the merger time, the galactic heights of 
the latter are expected to be considerably larger. 

Even in the absence of kicks (ns, ns) binaries, which 
escape mainly from dwarf galaxies, may form a common 
intcrgalactic field of these objects (Tutukov & Yungelson 
1994). It is suggested that most of star formation in the 
Universe occurs at low redshift (z <~ 1, see e.g. Madau 
1997, Connolly et al. 1997) and the main sites of star 
formation could be dwarf star-burst galaxies (e.g. Babul 
& Ferguson 1996). However, due to the large distances 
traveled by the (ns, ns) before they merge, a significant 
fraction merge at a large distance form the host galaxy in 
which the binary was formed. If 7-ray bursts are related 
to (ns, ns) mergers, a considerable fraction may occur at 
large distance from their host galaxy and they are not 
likely to be associated. 

4-3. Black hole plus neutron star binaries 

Our models suggest (Table 1) that the number of detected 
radio pulsars has to be at least 3 to 4 times larger than 
at present before one may hope to find a young pulsar ac- 
companied by a black hole. The only exception is in model 
H in the version in which hypercritical accretion onto the 
neutron star does not always result in swallowing of the 
companion. In this model one may expect the presence of 
2 or 3 (bh, ns) systems among the known young pulsars 
(see Table 1). The number of (bh, ns) binaries has then to 
be comparable to the number of HMBPs. Both these con- 



clusions are in conflict with the observations. The merger 
rate of (bh, ns) binaries is typically more than an order 
of magnitude smaller than that of (ns, ns) binaries. An 
exception to this is model D with small ?7coA, which pre- 
vents the formation of (ns, ns) systems by merging their 
progenitors in the common envelope phase. 

For gravitational waves detectors the signal-to-noise 
ratio is (Bonazzola & Marck 1994) 

' (Mr +M 2 ) 1 /6 v ; 

where are the masses of coalescing objects and D 

is the distance to them. E.g., for the range of black hole 
masses 4 to 17 M Q suggested by black-hole candidate bina- 
ries, one would expect that for (bh, ns) mergers detection 
distance may be 1.5 to 2.5 larger than for (ns, ns) mergers. 
Hence, the rates of detectable (ns, ns) and (bh, ns) merg- 
ers are expected to be comparable (about 1 / yr) . This has 
to be considered in the context of expected coalescence 
waveforms from the first detections (Thorne 1996). Also, 
the characteristics of 7-ray bursts expected from (bh, ns) 
mergers may differ from the ones from (ns, ns) mergers 
(Paczyhski 1991). 

Due to the large orbital period of (bh, bh) binaries, 
their expected merger rate is negligible and this number is 
not given in Table 1. A rate of detectable (bh, bh) mergers 
comparable to the rate of (ns, ns) mergers was predicted 
by e.g. Tutukov & Yungelson (1993a) and Lipunov et al. 
(1997). According to our current understanding this is a 
result of underestimating the effect of high mass loss in 
the stellar winds of stars more massive than ~ 40 M . 

5. Conclusions 

We discussed in detail the formation and evolution of bi- 
naries which contain two neutron stars or black holes. 
The orbital parameters of present-day population in the 
Galaxy are studied. For a successful model we require that 
also the results of the preceding evolutionary stages (oc- 
currence rates for supernovae and the formation of sin- 
gle pulsars, number of Be and high-mass X-ray binaries) 
comply with observations. As an additional restriction we 
require the observed high-mass binary pulsars to fit into 
the modeled P or b ~ e diagram for the appropriate age of 
the population. We summarize our conclusions: 

1. Deposition of orbital energy into common envelopes 
has to be very efficient for in-spiraling neutron stars, 
and sources other than orbital energy have to be in- 
voked. 

2. A kick velocity distribution with a high velocity 
tail and the bulk in small velocities (as given by 
Paczyhski 1990 with the parameters from Hansen & 
Phinncy 1996) gives a model which is in agreement 
with the observations. The same velocity distribution 
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provides a satisfactory model for single pulsars (Hart- 
man 1997). Increasing the fraction of low velocities 
may provide an even better fit to the observations. 

3. The best fit between the observations and the models 
is at an age of the population between several 100 Myr 
and 1 Gyr. This suggests that the majority of the re- 
cycled neutron stars are born with low magnetic fields. 

4. In the model of the population of (ns, ns) bina- 
ries which we consider satisfactory, their birthrate is 
— 3.4 x 10~ 5 yr _1 (assuming 100% binarity). This is 
consistent with the birthrate derived from observed bi- 
nary pulsars. The merger rate of (ns, ns) binaries in 
this model is <~ 2 x 10~ 5 yr _1 . The merger rate of (bh, 
ns) binaries is smaller by at least an order of magni- 
tude, but they may contribute significantly to the rate 
of expected gravitational wave detections because of 
the higher total mass of the systems. 

5. Binary neutron stars merge typically within 1 to 2 Gyr 
after formation. Hence, elliptical galaxies and the old 
stellar population in spiral and irregular galaxies may 
not contribute significantly to the rate of gravitational 
waves detections. When the detectors become sensitive 
to the events at cosmological distances this population 
may give an important contribution. 

6. The rate of (ns, ns) mergers is consistent with the 
estimated rate of 7-ray bursts if the latter are beamed 
into an opening angle of ~ 10°. 

7. The scenario for the evolution of close binary stars 
which assumes hypercritical accretion and the subse- 
quent collapse of the neutron star into a black hole in 
the common envelope is able to reproduce binary neu- 
tron stars with large orbital periods (~ 10 days), but 
fails to reproduce the short period systems. 

8. In the scenario in which a neutron star accretes hyper- 
critically and the binary avoids merging, the Galactic 
population of binaries which contain a black hole and 
a neutron star might well exceed the number of double 
neutron stars. The merger rate of these binaries is then 
approximately 3 x 10 yr . 

9. The number of observed single radio pulsars has to 
increase to about 2000, before one can hope to find a 
young radio pulsar accompanied by a black hole. 

10. A considerable fraction of the binary neutron stars es- 
cape from the tidal field of the Galaxy and travel up 
to <~ Mpc distances before they coalesce. 

1 1 . Short period neutron star binaries are, due to the high 
velocities received upon the supernova, virtually all 
ejected from the shallow potential well of globular star 
clusters. High-mass binary pulsars in globular clusters 
are therefore likely to be formed by a dynamical inter- 
action. 
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Fig. 3. The probability distribution for the number of Galactic 
(ns, ns) binaries younger than 100 Myr. The kicks are selected 
from a Maxwellian velocity distribution (model C in Table 1). 
The gray scaling is given in the upper left corner and is chosen 
to be the same as for the upper panel in Fig. 2. 
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Fig. 2. The probability distribution in model B for the present 
day orbital parameters of the Galactic disc neutron star bina- 
ries younger than 100 Myr (upper panel), 1 Gyr (middle panel) 
and 10 Gyr (lower panel). The gray scaling in this graph rep- 
resents numbers in the Galaxy. The darkest shades in the up- 
per, middle, and lower panels correspond to a total of Galactic 
23, 110, and 1100 (ns, ns) binaries with given combination of 
P or b and e. The "•" symbols in the upper panel represent the 
observed (ns, ns) binaries, the "★" symbols mark the position 
of the a system which might originate from a dynamical en- 
counter (the left *) and a system of which one companion is 
possibly not a neutron star (see Table 2). 
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Fig. 4. The probability distribution for model D (ri ca \ = 0.5, 
and kick according to Eq. 5) at the age of the population 
younger than 100 Myr. The gray scaling, given in the upper 
left corner, is chosen independently of the other models and 
gives 8.9 (ns, ns) binaries in the Galaxy for the darkest shade. 
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Fig. 5. The probability distribution for the Galactic number 
of (ns, ns) binaries younger than 100 Myr (upper panel) and 
10 Gyr (lower panel) in model H. The normalization gives 6.3 
(ns, ns) for the darkest shades in the upper panel and 522 for 
the lower panel 
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Fig. 6. The cumulative distribution of the center-of-mass ve- 
locities for neutron star binaries which merge within 10 Gyr 
from models A (dash-dotted line to the right), B (solid line) 
and C (dashed line). For model A, for comparison, the line to 
the left represents the velocity distribution of all (ns, ns) bina- 
ries and includes a rather large population of systems which do 
not merge within 10 Gyr. The dotted line represents the input 
kick- velocity distribution for model B (see Eq. 5). 
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Fig. 7. Upper panel - Normalized distribution of (ns, ns) bina- 
ries over merger ages in model B. The limit of 10 Gyr is adopted 
as the age of the galactic disc. Lower panel - The cumulative 
distribution of the (ns, ns) binaries over merger ages for the 
models A (dash dotted line), B (solid line) and C (dashed line). 
The fraction of mergers reached within 10 10 yrs indicates that 
a considerable fraction of systems does not merge within a 
Hubble time. 
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Fig. 8. The cumulative distribution of (ns, ns) binaries over 
the distance traveled before coalescence occurs. The line styles 
represent the results from the same model computations as 
denoted in Fig. 7. 
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Fig. 9. The cumulative distribution of (ns, ns) binaries over 
the galactic height z reached before coalescence occurs. The 
galactocentric distance at which the merger occurs is computed 
using Eq. 7. The line styles are as in Fig. 6. The subset of (ns, 
ns) binaries from model A which merge within 10 Gyr (right 
dash-dotted line) reach higher z- values then those with orbital 
period so large that they do not merge within the age of the 
Universe (left dash-dotted line) . 



